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Solving crystal structures from electron diffraction patterns rather than X-ray diffraction data is 
hampered by multiple scattering of the fast electrons within even very thin samples and the difficulty 
of obtaining diffraction data at a resolution high enough for applying direct phasing methods. This 
letter presents a method by which the effect of multiple scattering is being used for solving the 
phase problem, allowing the retrieval of electron structure factors from diffraction patterns recorded 
with varying angle of incidence without any assumption about the scattering potential itself. In 
particular, the resolution in the diffraction data does not need to be sufficient to resolve atoms, 
making this method particularly interesting for electron crystallography of 2-dimensional protein 
crystals and other beam-sensitive complex structures. 
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Electron crystallography becomes the method of choice 
when trying to determine the atomic structure of small 
crystalline volumes for mainly 2 reasons: a) electrons can 
be focused into sub-nanometer probes, and b) their scat- 
tering strength is about 4 orders of magnitude greater 
than that of X-rays. For structure determination of in- 
organic materials the strong interaction of electrons with 
the scattering medium is traditionally being viewed as a 
nuisance because existing methods for solving the phase 
problem in X-ray crystallography assume kinematic scat- 
tering, making these methods generally not directly ap- 
plicable to electron diffraction intensities. Solving the 
structure of 2-dimensional protein crystals from electron 
diffraction data alone, on the other hand, is hampered by 
that fact that these beam-sensitive samples do not gen- 
erally provide data of sufficient resolution to satisfy the 
preconditions for kinematic phasing methods. 

Recently, approaches which use the sensitivity of mul- 
tiply scattered X-rays to structure factor phases have 
proven useful P, 0, bl . The comparatively flat Ewald 
sphere describing the scattering of high energy electrons 
usually causes more than just 3 reflections to be close 
to the Bragg condition, especially in complex structures, 
making these methods generally inapplicable for phase 
determination from electron diffraction data. Also in the 
field of low energy electron diffraction (LEED) dynamical 
scattering is being applied for solving the phase problem 
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Without loss of generality it is assumed that the spec- 
imen surface normal is parallel to the zone axis as well 
as the z-axis, as illustrated in figure [1] The deviation 
of the 3-dimcnsional wave vector of the incident electron 
beam k = (k x , k y , k z ) = (k t ,k z ) from the z-axis is then 
usually quite small. Based on the Bloch wave solution 
to the Klein-Gordon equation of high energy scattering 
of a fast incident beam electron by the crystal potential 
[H, Q the intensity I gn {kt) of any point in the diffraction 
pattern defined by the reciprocal lattice vector g n of the 
particular diffraction spot and the tangential component 
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FIG. 1: Diagram illustrating the relevant parameters in real 
(a) and reciprocal space (b). a) In convergent beam elec- 
tron diffraction (CBED) experiments diffraction patterns for 
all incident beam directions whose fc t -component lie within 
a disc are recorded, b) The reciprocal space representation 
illustrates how changing the incident beam direction causes 
the Ewald sphere (dashed arch) to intersect the crystal poten- 
tial reciprocal lattice rods (REL rods) with varying excitation 
errors. 



of the incident electron wave vector k t is given by the 
modulus squared of the element of the scattering matrix 
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in the n row and the "central beam" column m, i.e. 
the beam index m for which the reciprocal lattice vec- 
tor g m = 0. Here A is a square matrix with the po- 
tential energy terms A n ^ m = Ug n _g m as its off-diagonal 
elements and the (relativistically corrected) terms = 
— {\g n \ 2 + 2<?„ ■k)/ r y related to the kinetic energy part of 
the modified Schroedinger equation along its diagonal. 

The scalar parameter T = irj\t depends on the wave- 
length A = (Aq~ 2 + Uo)~ 1 ^ 2 ~ A of the incident electron 
beam corrected by the mean potential Uq of the crystal, 
the specimen thickness t, and the relativistic correction 
factor 7 = 1 + \e\v/moc 2 (here mo and e are rest mass and 
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charge of the electron and v is the accelerating voltage). 

Note that I did not multiply the electron struc- 
ture factors Uff n -g m — 2m |e|V^ n _^ m //i 2 (Vg = 
il^ 1 J cell V{r) exp(2nig ■ r)d 3 f are the Fourier coefficients 
of the crystal potential V(r), where the integration is over 
one unit cell, and Q is the unit cell volume) in A with the 
usual relativistic correction factor 7 but included it in T 
instead, in order to separate true material constants (JJ$) 
from variable experimental parameters (A(u), j{v), t, kt). 
The additive constant of Uq for every element along the 
diagonal of A has also been omitted because it only pro- 
duces a general attenuation of the diffraction pattern due 
to its imaginary part and a global phase offset due to 
the real part of Uq, both immeasurable in a conventional 
diffraction experiment. 

The crystal potential V(r) — V r (r) + iF l (r) consists 
of a real part describing the elastic scattering and an 
imaginary part accounting for inelastic scattering pro- 
cesses which produce a non-isotropic attenuation in the 
scattered signal, if we collect the zero-loss scattered elec- 
trons only, as is routinely done in today's quantitative 
(energy filtered) electron diffraction experiments. This 
means that even for centro-symmetric crystals the Ug 
are in general complex, and have a non-zero phase. 

It has been shown that the Bloch wave expression (JTJ 
can be expanded in a "scattering path" or Born series 
[Icj |. providing a very general expression which agrees 
with those derived earlier by analyzing the multiple scat- 
tering process directly [HI, [H, Reference [1(3] pro- 
vides a recursive, but finite explicit expression for the 
coefficients ti t m (T, k t ) in the expansion 
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for arbitrary scattering path lengths, even in the (degen- 
erate) case of multiple excitation of the same reflection. 
The interested reader is referred to the original publica- 
tion where the complete expression along with its 
derivation is provided. 

The complex product of any two scattering 
path coefficients c n,lT,—,l gi -i,hi,—,h 92 -i 
Cl\,.., qi _ um {TM) ■ C^ hi _ hq2 _ um (T,k t y de- 
fines scale and phase factor with which any 
51 + 52 = q-monomial of structure factors 
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to the diffraction intensity in the beam g n . 

The convergence of above scattering path expansion 
d3|) depends on the value of T = wyXt and the modu- 
lus of the largest structure factor. Increasing the speci- 
men thickness or decreasing the accelerating voltage will 



therefore increase the largest order of significant monomi- 
als in the expansion. We will let q max be the length of the 
longest scattering path - or, in other words, the largest 
order of monomial - included in the expansion from here 
on. The coefficients c qi ' q2 = c qi ' q2 (T, k t ) (qi+q 2 < q max , 
I will omit the functional dependence for increased read- 
ability from here on) depend in addition to the speci- 
men thickness and accelerating voltage also on the easily 
variable experimental parameters kt- Measuring diffrac- 
tion intensities for a set of different values of incident 
electron beam tilt kt-, as in a convergent beam electron 
diffraction (CBED) experiment, allows a system of cou- 
pled multi-variate polynomial equations of degree q m ax 
to be defined whose solution is the set of structure fac- 
tors which best (in a least squares sense) describe the 
observed intensities for the set of experimental parame- 
ters. The resulting system of equations can be solved by 
the reformulation-linearization technique 14] which has 
been proven to converges to a global optimum. 

In order to keep the size of the set of polynomial equa- 
tions within manageable limits T should be small enough 
to let the first term in the expansion, the kinematic scat- 
tering intensity at reflection g n , m, 
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(S 3n = |fn/7 is also called the excitation error) be 
the dominant contribution, at least to the strong reflec- 
tions. Although the kinematic scattering approximation 
j( qmax — 2 ) is not sensitive to any structure factor phase 
invariants, it allows us to define sensible limits on the 
structure factor amplitude. 

It also provides a very simple way of determining the 
specimen thickness, the only unknown parameter in the 
c <?i,92 coefficients (accelerating voltage and incident beam 
tilt can usually be controlled very precisely). By fit- 
ting the (normalized) experimental 2-dimensional rock- 
ing curves Ig n (k t ) for all reflections with the kinematical 
( 2) — * 

one IgJ(T,k t ) given by expression (J3J), both, the speci- 
men thickness as well as a first estimate for the structure 
factor amplitude can be obtained. The form of expression 
([3]) allows this fit to be implemented very efficiently as a 
simple one-parameter search in the specimen thickness t. 
The structure factor amplitudes are then determined di- 
rectly for each trial value of t. This is, by the way a much 
more precise measurement of structure factor amplitudes 
than using the integrated rocking curve as is done in pre- 
cession electron diffraction [15| . 

Setting q max = 3 the scattering intensity at reflection 
9n (n 7^ to) becomes 
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The monomials C/g„-g, f/J can be written as the prod- 
uct of the 3 structure factor amplitudes and a phase fac- 
tor whose argument is the 3-phase invariant 4>g n -g t +4>g t — 
4>g n , indicating that IgJ(T,kt) must also be sensitive to 
these phase invariants. 

The reformulation-linearization technique developed 
by Shcrali and Tuncbilek [14j converts the non-linear 
polynomial equations into linear ones by treating each 
monomial as an independent variable. By defining ad- 
ditional linear equations that result from obvious rela- 
tionships between the monomials as well as inequalities 
produced by upper and lower bounds on the variables 
the polynomial set of equations is converted into a set of 
linear equalities and inequalities for which a solution can 
be found, if one exists, using standard linear program- 
ming algorithms. Once a solution has been found, the 
limits on the variables may be readjusted ("branch and 
bound" ) , resulting in a new linear programming problem 
and eventually the global solution. 

If the initial estimate of the structure factor ampli- 
tudes is fairly reliable, the limits within which the correct 
structure factor amplitude must lie may be defined very 
tightly. It is then very likely that the correct solution 
will be found already in the first linearization step, as 
has been the case for the test case described below. 

In order to verify the performance of the phase retrieval 
method described above it has been tested using simu- 
lated electron diffraction data. Since it has already been 
shown that Friedel asymmetries in dynamical diffrac- 
tion patterns of non-centrosymmetric structures can be 
used to solve the phase problem [3], I used simulated 
diffraction patterns of Si in the (110) projection, a simple 
centro-symmetric structure for demonstration purposes. 
If absorption is neglected, as has been done for this test 
case, kinematic and higher order scattering contributions 
to the diffraction intensities cannot be separated simply 
on the basis of their symmetry, as it is partly the case 
for non-ccntro-symmetric structures flit ] , making the re- 
construction of centro-symmetric structures even more 
challenging. 

Recently, a TEM-controlling software has been devel- 
oped in our laboratory (and tested on a Zeiss EM 912 and 
the Zeiss SESAM [18|) that allows automated acquisi- 
tion of energy-filtered pseudo-large angle CBED (pseudo- 
LACBED) patterns to be recorded on either slow-scan 



CCD camera or imaging plates |17[. The advantage 
of pseudo-LACBED patterns over conventional CBED 
patterns is that, although the radius of each individual 
diffraction disc may be several reciprocal lattice vectors, 
discs produced by the different diffraction spots do not 
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FIG. 2: a) Simulated LACBED disc for g = [113]. b) The 
coefficient e£|™(t = 3.985nm, k t ) given by expression [3] (un- 
sealed kinematic rocking curve). The square root of the pro- 
portionality constant between the 2 plots shown in a) and 
b) provides a first estimate of the structure factor amplitude 
for this reflection, c) Simulated LACBED disc for the kine- 
matically forbidden reflection [222]. Comparing the diffrac- 
tion intensity of a forbidden reflections with the correspond- 
ing kinematic rocking curve (d) allows its identification to be 
automated. 



overlap and may therefore be recorded in a single ex- 
posure. This apparently unphysical effect is produced by 
scanning (and, below the specimen, partially descanning) 
the angle of incidence of a nearly parallel beam during 
the exposure of the recording medium. The large disc ra- 
dius is necessary to observe fluctuations in the diffraction 
intensities across the discs despite the low sample thick- 
ness required to keep dynamical diffraction effects at a 
moderate level. For light matter (e.g. protein crystals) 
the sample thickness may actually be quite high, dras- 
tically reducing the need to use very large convergence 
angles. 

LACBED disc intensities for the following 27 beams 
have been computed by Bloch wave simulation: [222], 
[222], [002], [331], [331], [224], [224], [113], [113], [004], 
[220], [111], [111], [000], [111], [111], [220], [004], [113], 
[113], [224], [224], [331], [331], [002], [222], [222]. An 
accelerating voltage of 200 kV, a specimen thickness of 4 
nm and a disc radius of 3° have been used, a beam tilt 
angle which is easily achievable in actual experiments 
with the acquisition method described in the previous 
paragraph. 

Examples of the resulting diffraction data are shown 
in figure [5] a and c. While the high symmetry of centro- 
symmetric structures removes potentially useful Friedel 
asymmetries in the diffraction data, it also imposes sym- 
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a) Original V(r) b) Reconstruction c) Reconstruction 
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FIG. 3: Reconstructed real-space potential distribution 
within a single unit cell of Si in the (110) projection, a) Poten- 
tial distribution used to simulate the diffraction data, b) Re- 
construction from a scattering path expansion with q ma x = 3, 
i.e. one scattering order beyond the kinematic approximation. 
It is quite obvious that this order of approximation is not suf- 
ficient, c) Reconstruction from a scattering path expansion 
with q m ax = 3 showing excellent agreement with the original. 
Since the value of Uo cannot be reconstructed from the avail- 
able data, the potential is displayed with its minimum value 
set to in all 3 figures. 



metry constraints on the structure factors themselves. 
The 26 structure factors corresponding to the reflections 
listed above are reduced to 8 independent real-values 
variables (Uo was set to fro reasons mentioned earlier). 

In step 1 of the reconstruction the specimen thickness 
and initial estimates of the structure factor amplitudes 
were obtained by fitting kinematic rocking curves (see fig- 
ure[5]b and d) to the simulated pseudo-LACBED discs by 
nonlinear single parameter least-squares fitting w.r.t. the 
specimen thickness as described earlier, automatically re- 
moving forbidden reflections from the fit based on their 
poor agreement with the shape of the kinematic rocking 
curves. The best matching thickness was found to be 
3.985 nm, deviating only 0.3 from the thickness used for 
the simulation, most likely due to dynamical effects. 

In step 2 of the reconstruction the previously deter- 
mined sample thickness of 3.985 nm was used to define 
the c q ^ q2 (T, k t ) coefficients. Both, the definition of these 
coefficients as well as the setup of the linearized system 
of polynomial equations is fully automated, requiring as 
only input the specimen thickness, the reciprocal lattice 
vectors to be included in the reconstruction, the devia- 
tion from zone axis k t for each point within the LACBED 
discs and the diffraction intensity at that point. All of 
these parameters can be obtained directly from the (cal- 
ibrated) diffraction pattern. 

In step 3 of the reconstruction the previously defined 
system of linear equations and inequalities was solved 
using the lsqlin function provided by the Matlab opti- 
mization toolbox, which applies an active set method for 
solving linear programs in its medium-scale version. Re- 
peating the reconstruction for q rnax = 3 and q rnax = 4 
resulted in the projected potential plots shown in figures 



[3}3 and[3j: respectively. The reconstruction for q max — 4 
involved 33 different monomials, 11670 equalities and 122 
inequalities and took about 1 second on a 2.2 GHz lap- 
top equipped with 500 MB RAM to compute, solving for 
magnitude and sign of 6 distinct structure factors (being 
equal to zero, the structure factors of the forbidden re- 
flections had been removed automatically from the fit in 
step 1 of the reconstruction). 

The agreement of the second reconstruction including 
scattering paths up to length 4 with the original potential 
is excellent, demonstrating that this method works very 
well, provided that scattering paths of sufficiently high 
order are included. Increasing the number of reflections 
increases the number of unknowns, but also the number 
of available data. However, the number of monomials, 
equalities and inequalities increases as the q^ ax power 
of the number of structure factors to be determined. It 
will therefore be necessary to define a scheme by which 
the dominant monomials and the linear relationships be- 
tween them will be kept, discarding those which are not 
important - a path currently being pursued. 

In summary, a new method for solving the inversion 
problem of dynamical scattering from (convergent beam 
electron) diffraction data from moderately thin samples 
has been presented. Apart from the usually well-defined 
accelerating voltage the data required by the reconstruc- 
tion may be directly obtained from the diffraction pat- 
tern. No assumptions are being made about the scatter- 
ing potential within a single unit cell, making this method 
much more general than direct phasing methods applied 
to kinematical diffraction data. The method uses dy- 
namical scattering effects in the diffraction data to solve 
directly for the structure factor amplitudes and phases. 

Although the application of this method to X-ray 
diffraction patterns is quite conceivable, it is expected 
that it's main application will be in the field of structural 
electron crystallography. Simulated pseudo-LACBED 
data of Si (110) has been used to test this method. 

I want to thank Prof. H. D. Sherali for his very helpful 
suggestions on how to solve multivariate polynomial sets 
of equations. 
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